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ABSTRACT 

There are mainly two different approaches to measure the cosmic star formation his- 
tory: direct star formation rate density (SFRD) and steUar mass density as functions 
of redshift. Compilations of current observations seem to show a disparity in the two 
quantities, in the sense that the integral of SFRD is higher than the observed p^, (af- 
ter considering gas recycling). Using cosmological smoothed particle hydrodynamics 
simulations based on the concordance A cold dark matter model, we show that the 
two quantities become more consistent with each other when we consider the observed 
galaxy mass limit. The comparison between simulations and (dust corrected) observed 
cosmic SFRD shows a good agreement, while the observed is significantly lower than 
the simulation results. This can be reconciled if the current high-z galaxy surveys are 
missing faint low-mass galaxies due to their flux limit. Our simulated GSMFs have 
steep low-mass end slopes of a < —2 at z > 3, and when these numerous low-mass 
galaxies are included, the total p» matches with the integral of SFRD. 

Key words: methods: N-body simulations — galaxies: evolution — galaxies: forma- 
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1 INTRODUCTION 

The cosmic star formation (SF) history is a fundamental 
quantity that illustrates the buildup of galaxies. It can be 
estimated by direct measurement of SFRD as a function of 
redshift through rest-frame Ha, UV, and far-IR emission 
with dust extinction corrections (e.g., iHopkins fc BeacomI 
I2OO6I ). The stellar mass density (p,) at various redshifts is 
an indirect measurement of the cosmic SF history, because 
it is indeed the integral of SFRD, after considering the recy- 
cling of gas into interstellar medium. Given this simple rela- 
tionship, in principle the measurements of the two quantities 
should be consistent with each other. However, in this paper 
we show that this is not the case for current observational 
estim ates, and other authors have reached similar conclu- 
sions (^Nagamine et al."2004VOuchi et a"Lll2004l : Ivan DokkumI 
[2OO8. : ,Wilkins ct al. 2008; Dave 2008). 

Over the past decade, large samples of color selected 
high-z galaxies such as Lyman break galaxies (LBGs) have 
enabled us to measure the cosmic SFRD at 3 < z < 6 
from rest-frame UV luminosity density (e.g., ISteidel et al.l 
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119991 ). A number of near-infrared observations have also 
constrained the galaxy stellar mas s function (GSMF) at 
z <4 (e.g.. iMarchesini et alllioogl ). Furthermore, the ad- 
vent of the Wide Field Camera 3 {WFC3) on board the 
Hubble Space Telescope (HST) has dramatically improved 
our ability to measure the rest-frame UV light from galaxies 
at z > 6. Combined with the measurement of rest-frame op- 
tical light by the Infrared Array Camera {IRAC) on Spitzer 
Space Telescope, it is now possible to measure the charac- 
teri stic stellar mass of galaxies at z > 6 using stacked spec- 



tra jYan et al.ll2010l:lBouwens et ai]|2010l : iLabbe et all2010l : 



iGonzalez et al.ll2010l V 

However, the sources at z > 6 are very faint, and so 
far we have only detected the massive end of GSMF. In 
addition, there are still significant uncertainties in the esti- 
mates of SFRD and p*. Given this situation, it would be 
useful to obtain predictions on SFRD and p* from the- 
oretical models. In particular, cosmological hydrodynamic 
simulations have b een widely used to investigate cosmic 
star formation ( e.g..lCen fc Ostrikeijll992l: iKatz et al.lll996l: 



Springcl fc Hernguistl l2003bl : iNagamine et al.l I2OO6I : 



2008i : ,Schave et al.ll20ld ) 



Davel 



In this Letter, we focus on the cosmic SF history at 
z > 2 using cosmological simulations. The remaining of this 
paper is organized as follows. In Section [2l we describe the 
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Name 


Box-size 


AT 




^gas 


e 




N216L10 


10.0 


2 X 216^ 


5.96 X 10"^ 


1.21 X 10^ 


1.85 


2.75 


N400L34 


33.75 


2 X 400^ 


3.49 X 10^ 


7.31 X 10^ 


3.38 


1.0 


N600L100 


100.0 


2 X 600^ 


2.70 X 10* 


5.66 X lO'^ 


4.30 


0.0 



Table 1. The three different resolution and volume simulations employed in this Letter. The box-size is given in units of h~^Mpc, A'^p 
is the particle number of dark matter and initial gas (hence X 2), oidm and mgas are the masses of dark matter and initial gas particles 
in units of h~^MQ, respectively, e is the comoving gravitational softening length in units of /i~^kpc, and z^^^ is the ending redshift of 
the simulation. Note that the mass of star particle is a half of the initial gas particle. The value of e is a measure of spatial resolution. 
The name of each simulation is based on the particle count and its box size. 



cosmological hydrodynamic simulations. In Section [S] we 
construct the composite GSMF by combining the samples 
of galaxies from simulations with different resolution and 
volumes. We study the cosmic SFR and p* evolution in Sec- 
tion [l] We summarise and discuss our findings in Section [S] 



conversion is necessary to obtain SFRD from our simula- 
tions. This is why our simulations are helpful to resolve this 
problem between SFRD and p*. 



2 NUMERICAL SIMULATIONS 

We use the modified version of the tree-particle-mesh 
smoothed particle hydrody namics (SPH) c ode GADGET- 

3 (originally described in ISpringej l2005h . Our conven- 
tion al code in cludes radiative cooling by H, He, and met- 
als l|Choi fc Na gamine 20 09i), heating by a unifor m UV 
backgr ound of a modifi ed Haardt fc Madaul l|l996l) spec- 
trum l|Katz et all Il996l : iDave et al.l Il999l ). star forma- 
tion, its feedback, a phenomenological model for galac- 
tic winds, and a s ub-resolution model of multiphase ISM 
l|Springel fc Herriguist 2003a ). For the star formation model, 
we use the "Pressure model" which reduces the high-z SF RD 
l|Schave fc Dalla Vecchiall2008l : IChoi fc Nagamine''2 010l) rel- 
ative to the previous model by Springol fc H ernguistl 
iioosi). For the galactic outflow model, we use the Mul- 
ticomponent Variable V elocity wind model developed by 
IChoi fc Nagamin^ l|201ll) . We adopt the following cosmolog- 
ical par ameters that are cons istent with the WMAP best-fit 
values l|Komatsu et al.1 120111 ): fi™ = 0.26, = 0.74, Qt = 
0.044, h = 0.72, ns = 0.96, and as = 0.80. To identify galax- 
ies in simul ations, we use simplified variant of the SUBFI ND 
algorithm jSpringel et al ] |200ll : IChoi fc Nagamine! l2009l ) . 

Although cosmological simulations have been widely 
used to study galaxy formation, there are inevitable lim- 
itations due to resolution and box size. One critical limi- 
tation is that the low-mass galaxies are not captured in a 
low-resolution simulation with a large box size, and the high 
mass galaxies are missed in a small box size simulation due 
to limited box size. To alleviate this problem, we construct 
composite GSMFs in a wide range of galaxy stellar mass us- 
ing a few simulations with different resolution and volumes 
(see Table [Ij. 

An important point to note is that the SFRD is one of 
the most direct output of our ab initio cosmological simu- 
lations; our simulation follows the gas and dark matter dy- 
namics within the framework of ACDM model, and converts 
part of the gas into star particles when the SF threshold den- 
sity is satisfied at every time step. Therefore the SFRD is 
the most direct outcome of gas dynamics, and no uncertain 



3 GALAXY STELLAR MASS FUNCTION 

Figure [T] shows the composite GSMFs for different redshifts. 
Three simulations with different resolutions and volumes 
(N216L10, N400L34, and N600L100 runs) are used to con- 
struct the composite GSMFs, which cover the galaxy stellar 
mass range of lO^ft'^M© < < lO^^/i'^Mg. The low- 
est mass galaxies are represented by the N216L10 run, the 
intermediate mass galaxies by the N400L34 run, and the 
most massive galaxies by the N600L100 run. The covered 
mass range for each run is determined by the following two 
criteria: First, the lowest M* of a given simulation is equiv- 
alent to that of 32 star particles, which is ^ 10^/i~^Mq, 
~ lO^h'^AlQ, and ~ lO^fe-^M© for the N216L10, N400L34, 
and N600L100 runs, respectively. Second, ff the GSMFs 
overlap between different runs, the one from a larger vol- 
ume run is used, because a smaller volume run tend to un- 
derestimate the number density of galaxies at the massive 
end owing to the limited box size. Figure [1] shows marginal 
overlaps between two simulations for presentation purposes. 
Three segments of GSMF from three simulations smoothly 
form one mass function. 

The simulated GSMF shows a reasonable agreement 
with observations (yellow and blue shade) at 2 < 2 < 4, 
although the simulated GSMF is slightly lower than the ob- 
servations at the very massive end of M* > 10^^ Mq. 
The number density of these high-mass galaxies is so low 
that there are only about 1 — 10 galaxies in a volume of 
(100/i~^Mpc)^, therefore this discrepancy could be due to 
the cosmic variance. Another possible reason is that the 
lower mass galaxies are scattered up into the more mas- 
sive bins due to observational errors, causing it to be higher 
than the intrinsic simulated GSMF. 

A marked feature in the simulated GSMF is the very 
steep slope at the low-mass end. When fitted with a 
Schechter function, the simulated composite GSMF has a 
faint-end slope of a —2 at z = 3 for the mass range 
of 10^ h'^MQ < M* < W^°h~^MQ, while the observations 
sugg est a ^ —1.6 based on the data at M* > lO^'^ Mq 
(e.g. iMarchesini et al.1 120091 '). The simulation suggests that 
a large fraction of stars is located in the galaxies with A'h < 
10^° h'^MQ at z = 3. We find that 22% of total p* is located 
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Figure 1. Simulated GSMFs at 2=3, 4, 5, 6, 7, and 8, from our three simulations. The low mass range (M* < 1 X 10^h~^ AIq) is 
represented by the N216L10 run, the intermediate mass range by the N400L34 run, and the high mass range (M* > 1 X 10^h~^ Mq) 
by the N600L100 run. The shaded regions show the range of observations at 2 < z < 3 (blue shade) and at 3 < z < 4 (yellow shade) 
l|Marchesini et ahllioogl ). The black dotted lines show the low-mass end slopes of a indicated in each panel. The simulated GSMFs show 
good agreement with observational data, and they have significantly steep low-mass end slopes in the mass-scales where the observations 
are not available due to faintness of those sources. 



in the galaxies with > 10^°/i"^Mq, 34% in W^H'^Mq < 
< 10^°h-^MQ, 31% in IO^/i-^Mq < < IO^/i-^Mq, 
and 13% in lO^/i"^M0 < M* < 10**/i~U/q for the simu- 
lated GSMF ai z — 3. In general, current galaxy surveys 
cannot detect galaxies with < 10^'^ Mq due to their 
flux limit at z > 3, and deeper observations are needed to 
probe the GSMF at lower masses. 

Another notable feature of simulated GSMF is that its 
low-mass slope becomes steeper as z increases. For example, 
the slope is a ~ —2.4 at 2 = 6, and a ~ —2.9 at z = 8. We 
perform more rigorous fitting with a Sc hechter function in 
a separate publication (Ijaacks et al.ll201ll ). which still shows 
a < —2. This steepness of simulated GSMF has significant 
consequences when one measures the global quantities such 
as cosmic SFRD and p,. We note that the integration of 
GSMF to M* = will diverge if q < — 2, and it requires a 
low- mass limit to compute the total p*. In our simulations, 
the lowest galaxy mass is determined by the atomic cooling 
limit of T ~ 10* K, which corresponds to ~ lO'^M©. 
Therefore the summation of simulated galaxy stellar mass 
stops at M* ~ W'^h'^MQ. 



4 HIGH-Z SFRD AND p* 

Figure [2] shows the evolution of SFRD and p* as functions of 
redshift. We decompose the total SFRD and p* into different 
galaxy mass ranges using our simulations. We note that the 
galajcy masses used for decomposition are the values at each 
redshift. At z > 5.5, most of the stars form in low-mass 
galaxies with M* = lO'^ - lO'^ft-^M©, while at z < 3.5 most 
form in massive galaxies with A/* > Mq. The same 

trend is also found in the p, evolution. It shows that the low- 



mass galaxies are preferred sites for SF in the early Universe 
as one would expect in the hierarchical structure formation 
model. 

The black solid lines in Figure [2] show the sum of SFRD 
and p* in three different galaxy mass ranges; i.e., the total 
SFRD and p* for aU galaxies with M^, > 10''h~^MQ. An 
important point to notice in Figure [2] is that the simulated 
total SFRD is well within the observed range (yellow shade 
and the o ther data points inc luding those from gamma-ray 
bursts bv lKistler et al.l (|2009l ')'). while the simulated total p, 
is clearly higher than the observed p*. Since the black solid 
line in Fig. [2)d is a direct outcome of simulated SFRD in our 
simulation, the two are theoretically consistent with each 
other as we described at the end of Section O Whereas in 
the observations, different techniques are used to estimate 
SFRD and p* with different mass limits, therefore the in- 
consistency could arise between the two measurements. The 
comparison of black solid line in Fig. [2}d and the observed 
data clearly indicates the inconsistency between the obser- 
vational measurements of SFRD and p*. 

According to our simulation, the solution to the above 
discrepancy is that the current observations are missing the 
starlight from low-mass galaxies at 2; > 3, and that the low- 
mass end of GSMF is steeper than what current observa- 
tions suggest. This would explain why the observed p* is 
lower than the si mulation result. This po int was partly ad- 
dressed earlier bv lNagamine et al.l (|2004l ). but our new set 
of simulations reinforce their argument with a better handle 
on the GSMF over a wider range of galaxy mass. 

The uncertainty in the dust extinction correction needs 
to be addressed, as the rest-f r ame U V light can be extin- 
guished by dust. lOuchi et all (|2009l ') estimated the SFRD 
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Figure 2. Cosmic SFRD (left) and p» (right) as functions of redshift. In each panel, tlie contribution from three different stellar mass 
ranges (lO'^ — lO*ft.~^M0, 10® — 10^h~^MQ, and > 10^h~^MQ) is shown. The black solid line in each panel show the total cosmic 
SFRD and p<, from al l galaxies in three different mass ranges. Left: The yellow shade is the locus of the observed data compiled by 
iNagamine et al~ corrected for the Salpeter IMF. The black op en pentagons with error bars are the observational estimates of 

SFRD using the high-z from 7-ray bu rst data from lKistler et al ] l|2009h . The black open triangles with error bars are the compilation of 
observational estimates of SFRD from I Gonzalez et al. 1 1I2OIOI') for galaxies with > 0.3L*_^. The black solid pentagons with error bars are 
the compilation of observational estimates of SFRD from lBouwens et al.l l|2010l ) for galaxies with > 0.05L*_3 with dust correction. The 
black squares and in verted-triangl e s at z ~ 7 are the observational estimates of SFRD with (solid symbols) and without (open symbols) 
dust correction from [Ouchi et al] ||2009| ): the squares are the SFRD integrated down to L ~ O.IL* and the inverted- triangles are the 
SFRD integrated down L = 0, which agrees well with our total sum (solid black line). The black solid line shows good agreement with 
the observed range of SFRD shown by the yellow and other observational data points. Right: The yellow shad e is the estimated pt by 
integrating the GSMFs over t he stellar mass r ange of M* > 10®/i~^Mq (Marchesini ct al. 200^). The GSMFs in lMarchesini et al.l (12009 ) 
are established by fitting the Schechtcr (197d ) function to t he observed galaxies with A/*^ > 10^'^ Mq. The black open triangles with 
error bars are the observed data from Gonza lez et al. 1 ll2010h . The black solid triangle at z ~ 2 with error bars is the o bserved data from 
iReddv fc Stcidcl ( 200g) for L > 0.083L*_2 3. The black solid line is clearly higher than the observational estimate bv lMarchesini et al.l 
||2009| ). All observational data are corrected for the lSalpeterl lll955l) IMF. 



at 2: ~ 7 both with and without dust corrections as shown 
in the Figure [2^, integrating their estimated UV luminos- 
ity function down to zero luminosity. Their total SFRD es- 
timate (inverted solid triangle) agrees with our simulation 
very well. Therefore it is unlikely that the dust extinction 
correction plays a significant role in resolving the inconsis- 
tency between the observed SFRD and p* estimates. 

Furthermo re, the observational estimates by 

iGonzalez et al.l (|2010l ) at z > 4, which takes into account 
of galaxies with > 0.31/^^3 show a good agreement with 
our simulation results for galaxies with > 10^ Mq 
for both SFRD and p«. The correspondence between 
luminosity and stellar mass has a large scatter, but we 
find 1/2=3 ~ 10^" Mq both observationally and in our sim- 
ulatio ns. Therefore the > 0.3Z/J=3 limit of iGonzalez et al.l 
l|2010l) would correspond to roughly M^, ~ lO^ ' ^/t'^Ma. 
Observational estimates by iBouwens et al.l (|2010l ) include 
galaxies with an order of magnitude less luminous galaxies 
(> 0.05I/*=3, corresponding to roughly > 5 x 10** Mq), 
so it would correspond to the sum of a part of red line 
(10* < M*/Mq < 10^) and the blue line (M* > 10^ Mq). 
Therefore a proper comparison between simulations and 
observations with the same galaxy mass limit does not show 
the inconsistency of our concern, and the nice agreement 



bet ween simulation s and the data of lOuchi et al.l (|2009| ) 
and IGonzalez et al.l l|2010l ) support our argument that the 
missed low-mass galaxies can account for the discrepancy 
between the observed SFRD and p*. 

iJaacks et all l|201ll ) also demonstrated that the UV lu- 
minosity functions from our simulations and observations 
show an acceptable agreement, while the GSMFs show 
a larger discrepancy. Generally, SFRD is estimated from 
the rest-frame UV measurement, and the nice agreement 
between the SFRDs of simulation and observations sup- 
ports the agreement in the UV luminosity function as well. 
Whereas the p, is computed by integrating the GSMFs, 
which constrained less robustly at high-redshift due to un- 
certain light-to-mass conversion. This also supports our ar- 
gument that the origin of the discrepancy being in the low- 
mass end slope of the GSMF. 

The reasonable agreement between our total SFRD 
(black solid line in Fig. [2^) and the data points of 
iKistler et al.l l|2009l ) corroborates our argument that the low- 
mass galaxies are the dominant contributor to the total 
SFRD at high redshift. This is because the GRB obser- 
vations are independent of the rest-frame UV fiux limit of 
galaxy surveys, and GRBs would occur wherever the star 
formation takes place. However the GRBs may also be bi- 
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ased tracers of star formation, and the current GRB data 
suggests that an enhancement of the GRB rate at high-2 is 
necessary to account for the observed redshift distribution 
of GRBs jVirgih et al.ll201lh . A GRB rate increase due to 
lower metallicity at high-z may not be the singular cause, 
and other effect such as the evolution of the GRB luminos- 
ity function break with redshift might be important. Given 
these complications, further studies are needed in the future 
to better understand the relationship between the GRB rate 
and cosmic SFRD. 



5 DISCUSSION AND CONCLUSION 

Based on the results of self-consistent cosmological hydro- 
dynamic simulations, we argued that the current high-z ob- 
servations are missing low-mass galaxies with Af* < 10^ 
at z > 4 when accounting for the total p«, and this results 
in the inconsistency between the observational estimates of 
SFRD and p,. But are there any other possibilities to ex- 
plain the above inconsistency, or is it possible that our sim- 
ulations are incorrect? Below we discuss two points for such 
possibilities. 

First possibility is that our simulation might be over- 
producing the low-mass galaxies at z > 4. Recent observa- 
tional and theoretical studies show that the low-mass, high- 
z galaxies could have lower SF efficiencies than the nor- 
mal local galaxies (jWolfe fc Chenll2006l : I Gnedin fc Kravtsovl 
[2010 1 owing to lower molecular hydrogen fractions in low- 
metallicity environments. If this is true, our current simula- 
tions might be overpredicting the p« at z > 4. However if we 
revise our SF model to account for this effect, the SFRD will 
also decrease together with p*, and we might underpredict 
SFRD instead. In the future we will revise our SF model 
to consider the II2 fraction in high-2; galaxies, and evaluate 
how strong this effect would be. 

The second possibility is that the IMF at high- 2: could 
be differe nt from the loc al one. All data in Figure [5] as- 
sume the ISalpeteil (|l955h IMF. As we emphasized in the 
last paragraph of Section[2j the amount of gas converted into 
stars is a direct output of our simulations, therefore there 
are no uncertainties as to stellar IMF except the instanta- 
neous gas-recycling fraction, i.e., the ratio between the gas 
expelled from stars (via stellar winds and supernovae) and 
the total initial stellar mass. In our simulations, we assume 
the Salpeter IM F who se instantaneous gas-recycling frac- 
tion i s ~0.1 l|Spring cl & Hcrnquist 20033)- 

If we use the 

IChabried (I2003D IMF that has a lower number of low-mass 
stars, (3 increases to ~0.2, which would reduce the simu- 
lated p« by about 10%. Note that /3 here only considers the 
gas-recycling from massive stars. If we take into account of 
the contribution from long-lived, low-mass stars, the value 
of P will increase by a factor of a few. However, this pa- 
per focuses on the high-z galaxies, and we can safely ignore 
the gas-recycling from low-mass stars. Therefore the change 
from Salpeter to Chabrier IMF cannot fully account for the 
inconsistency between SFRD and p, . 

The top-heavy IMF in hig h-z galaxies has been spec- 
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ulated by many authors (e.g., iLarsonl 20051: Fardal et al.l 
: Ivan DokkumI l2008l : IWilkins et al.1 l20od : 
2010|), and it has a much higher value of /3. 



Consequently, the resulting p* will decrease significantly for 



a given SFRD, and it may alleviate the discrepancy in p*. 
However, the current theoretical models and observational 
suppor t for top-heavy IMF are not very robust; for ex- 
ample, Ivan Dokkum fc Conrovl l|2010l ) suggested abundant 
low-mass stars in high- 2: star-forming galaxies, which argues 
against a top-heavy IMF. 

In summary, we find a good agreement between our 
self-consistent cosmological simulations and the observa- 
tional estimates of SFRD and p, if we limit the compari- 
son to galaxies with > l(f Mq. In particular, the 
consistency b etween our simulati o ns and the o b servat ional 
estimates bv iKistler et al.1 (|2009l '). lOuchi et all l|2009l ') and 
iGonzalez et al.l (|2010l ')" at 2; > 4 is very encouraging. Our sim- 
ulations predict the existence of numerous low-mass galaxies 
with M* < 10^ Mq, and these low- mass galaxies are not 
included in the current observational estimates of p, . 

This paper demonstrates that the current observational 
estimates and understanding of the formation of high-z 
galaxies are still uncertain. To resolve this problem, we 
need future observations with increasing sensitivity and re- 
duced uncertainties (e.g., JWST) to provide more robust 
constraints for the abundance of low-mass galaxies at 2 > 4, 
as well as the improvement in the theoretical modeling of 
galaxy formation, particularly on star formation and its 
feedback. 
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